

setwd("/Users/italolopez/Dropbox/NEP Evaluation - The role of beliefs/JPE Final Version/JPE Submission/Replication package/Outputs/Tables")
############################################################################
############################################################################
### DV: zptevir_irtscore2 ###
############################################################################
############################################################################

################# MV: Behavioral Index ###################
set.seed(7688);
require(MASS)
a=0.108
b=0.318
acov <- matrix(c(
  0.002116, 0,
  0,0.0004
),2,2)
rep=20000
conf=95
pest=c(a,b)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = FALSE) 
M=mean(ab)
med=median(ab)

################# MV: Beliefs Index ###################
set.seed(7688);
require(MASS)
a=0.113
b=0.266
acov <- matrix(c(
  0.002209, 0,
  0,0.0004
),2,2)
rep=20000
conf=95
pest=c(a,b)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = TRUE) 
M=mean(ab)
med=median(ab)


################# MV: Behavioral and Beliefs Index ###################
set.seed(7688);
a1=0.108 #estimated coefficient a1
a2=0.113 #estimated coefficient a2
b1=0.245 #estimated coefficient b1
b2=0.164 #estimated coefficient b2
a1std=0.046 #SE of coefficient a1
a2std=0.047 #SE of coefficient a2
b1std=0.022 #SE of coefficient b1
b2std=0.021 #SE of coefficient b2
rep=20000 #number of simulated values
conf=95 #confidence level
a1vec=rnorm(rep)*a1std+a1 #create vector of simulated a1 coefficients
a2vec=rnorm(rep)*a2std+a2 #create vector of simulated a2 coefficients
b1vec=rnorm(rep)*b1std+b1 #create vector of simulated b1 coefficients
b2vec=rnorm(rep)*b2std+b2 #create vector of simulated b2 coefficients
total=(a1vec*b1vec)+(a2vec*b2vec) # simulated total indirect effects
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(total,low) #lower limit of confidence interval
UL=quantile(total,upp) #upper quantile for confidence interval
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = TRUE) 
M=mean(ab)
med=median(ab)

